% Process the eej strength derived from satellite data
% and compare it with ground deribed eej indices
%start date 30 OCT 2005
%latest date 29 MARCH 2006
clear;
more off;
load c:\manoj\projects\eej\OBS_array.mat OBS
%OBSdH1 description 
%The file OBS is Structured array with stations as listed below
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP,17 UJJ,18 KOD,19 YAP, 20 ]

% Note please update f:\manoj\projects\eej\used_obs.xls if more data are
% added to OBS_array.mat

%-------- correlate the CHAMP derived EEJ strength and ground delta H
% with LT correction
%--------------------------------------------------------------------

[a,b] = xlsread('c:\manoj\projects\eej\used_obs.xls');

bin5 = [-60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60];
bin10 = [-60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60]; 
bin20 = [-60 -40 -20 0 20 40 60];
obs_pairs = [2,3;4,1;5,6;15,14;19,20;18,3];
obs_title = ['ETT-HYB';'TIR-ABG';'HUA-FUQ';'AAE-ELT';'YAP-BIK';'KOD-HYB';]
    %Control Center
%[1-ABG,2-ETT,3-HYB,4-TIR,5-HUA,6-FUQ,7-MBO,8-GUI,9-AAE,10-QSB,11-GUA,
%12-CBI,13-PND,14-ELT,15-AAE(new),16-NGP,17-UJJ,18-KOD,19-YAP, 20-BIK ]
eej_data_set = 2;%1-old, 2-new
satdata  = 2; % %Pear Current Density = 1; av cur distan 2;
obslt    = [9.9 13.1]; %NOTE the LT are represented as for eg 8.99999 or 9.001 the filter
incrlat  = 10;
x1data   = bin10;
mindate  = tmjd(2000,01,01,0,0,0,1); % starting date
maxdate  = tmjd(2002,12,31,0,0,0,1)+1; % ending date
sqcor    = 0; %SQ correction
ltcor    = 1; %LT correction
robft    = 1; %Robust outlier elimination
plotopt  = 'b-'; %[b g r k]
plotcol  = 'b';
plotchar = 'b';
ncoeff   = 9;%same for hourly means 20 for minute means
kp_thrashold = 2;

switch eej_data_set,
    case 1,
        data = load('c:\manoj\projects\eej\eej_par_all.txt');
            switch satdata,
                case 1,
                eej_index = 9;
                case 2,
                eej_index = 10;
            end;
    case 2,
        load c:\manoj\projects\eej\EEJ_Stefan data;
        sqcor    = 0;
        switch satdata,
            case 1,
            eej_index = 5;
            case 2,
            eej_index = 6;
       end;
end;

for i_pair = 4:4,
x1data   = bin10;
eq_obs = obs_pairs(i_pair,1);
sq_obs = obs_pairs(i_pair,2);

%INFORMATION
%The EEJ parameters in the file
%1-mjd,2-longitude of eq-crossinfg, 3-zero of current position-south
%4 -north,5-half-width of current - south,6-north
%7minimum current south,8-north,9-pear current density 10-average
%difference 11 - diff to south minumum 12 - diff to north minum, 13-total
%east current ,14 total return current., 15 -f,16 - eq angle, 17-
%eq-direction, 18 - mgnteic equator latitude

    
%load f:\manoj\projects\eej\guacm4.mat index sq_est
        data1 = [];
        deltah = OBS(eq_obs).data - OBS(sq_obs).data;
        k = nanmean(reshape(deltah,[24,length(deltah)/24]),2)';  
%         k=circshift(k',[0,round(OBS(eq_obs).long/15)]); % 
        obs_index = 2;
        data1(:,1) = OBS(eq_obs).mjd;
        data1(:,2) = OBS(eq_obs).data - OBS(sq_obs).data;
%       load f:\manoj\projects\eej\ettcm4.mat index sq_est;
        xtime = 0.5:1:23.5;
% 		xtime = 1:24;
        obslon = OBS(eq_obs).long;
        if obslon > 180,
            obslon = -(360-obslon);
        end;
        st = obslon - (incrlat*6+incrlat/2);
        en = obslon + (incrlat*6+incrlat/2);
        
[p,s] = polyfit(xtime,k,ncoeff);%In case of minute data change the degree to 20
%----------------------------

ndd = zeros([1,25]);
nloop = 1;
lto = [];
ltc = [];
xdata = [];
nbindata = [];
ceff = [];
cerr = [];
nbin = 1;

for ki = st:incrlat:en,
L = data(:,2) > ki & data(:,2) < ki+incrlat & data(:,1) < maxdate  &  data(:,1) > mindate; %The last filter is to select only data before Dec 31, 2002
IndiaE = data(L,:);
ndata = 1;
ind1 = [];
ind2 = [];
ltt1 = [];
ltt2 = [];
nnodata=1;
for i = 1:length(IndiaE(:,1)),
    ind = findnearest(IndiaE(i,1),data1(:,1));
    [lt1,dummy1] = champ_lt(IndiaE(i,1),IndiaE(i,2),obslon);
    [dummy2,lt2] = champ_lt(data1(ind,1),IndiaE(i,2),obslon); % This is because observatory FDAY ~= CHAMP FDAY
    ltt1(i) = lt1;
    ltt2(i) = lt2;
    ltt3(i) = dummy1;
    ltt4(i) = dummy2;
    obst(i) = data1(ind,1);
    chpt(i) = IndiaE(i,1);
    
    if eej_data_set == 1,
        
         if lt2 >= obslt(1) & lt2 <= obslt(2) & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)),% & kp(ind) <= 2.5, %The last fliter take care the nodata (9999) points
            [ut_obs] = (data1(ind,1)-floor(data1(ind,1)))*24;
            [ut_sat] = (IndiaE(i,1)-floor(IndiaE(i,1)))*24;
             ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,ut_sat); %Getting the expexted EEJ strength from poly coefficients for champ 
            f_ob(ndata) = polyval(p,ut_obs); %Getting the expexted EEJ strength from poly coefficients for obs 
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
            ndata = ndata+1;
        else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
        end;
    elseif eej_data_set == 2,
            kp1 = IndiaE(i,10);
            cmax = IndiaE(i,5);
            if lt2 >= obslt(1) & lt2 <= obslt(2) & lt1 >= obslt(1) & lt1 <= obslt(2) ...
            & abs(data1(ind,obs_index)) < 300 & data1(ind,1) < max(data(:,1)) ...
            & kp1 <= kp_thrashold & cmax >= 0.03, %The second last fliter take care the nodata (9999) points
            [ut_obs] = (data1(ind,1)-floor(data1(ind,1)))*24;
            [ut_sat] = (IndiaE(i,1)-floor(IndiaE(i,1)))*24;
            ind1(ndata) = ind;
            f_ch(ndata) = polyval(p,ut_sat); %Getting the expexted EEJ strength from poly coefficients for champ LT
            f_ob(ndata) = polyval(p,ut_obs); %Getting the expexted EEJ strength from poly coefficients for obs LT
%             f_ch(ndata) = interp1(1:24,k,lt1,'spline'); %Found spline interpolation working better 18.1.06 !
%             f_ob(ndata) = interp1(1:24,k,lt2,'spline');
%             fprintf('Lon = %5.1f i = %d f_ch = %f f_ob = %f lt1 = %f lt2 = %f kp = %3.1f\n', ki,i,f_ch(ndata),f_ob(ndata),lt1,lt2,data1(ind,end));
             ndata = ndata+1;
             else,
            ind2(nnodata) = i;
            nnodata = nnodata + 1;
            end;
    end;
end;

if length(ind2) > 0,
    IndiaE(ind2,:) = [];
end;

if length(ind1) > 1,
    obsE = data1(ind1,:);
    X = IndiaE(:,eej_index);
    if ltcor == 1,
        Y = obsE(:,obs_index).*(f_ch./f_ob)';
    elseif ltcor == 0,
        Y = obsE(:,obs_index)';
    end;
    
    if robft == 1,    
%    Robust outlier rejection  -- START
     if length(X) > 3,
     [rs,rstat] = robustfit(X,Y);
     residuals = abs(rstat.resid);
     LLL = residuals == max(residuals); %Just remove the most noisy point
     X(LLL) = [];
     Y(LLL) = [];
     end;
    % Robust outlier rejection  --- END
    end;

    [cof,err] = corrcoef(X,Y);  % correlation & Correction for observatory data 
    ceff(nloop) = cof(2,1);
    cerr(nloop) = err(2,1);
    fprintf('%4.3f ',cof(2,1));
    xdata(nloop) = ki+incrlat/2;
    ndd(nloop) = ndata-1;
    lto(nloop) = min(ltt2);
    ltc(nloop) = max(ltt2);
    nloop = nloop +1;

    %----------- PLOT SCATTER
%subplot(2,3,nloop-1);
if x1data(nbin) == 0,
figure(2);
subplot(3,2,i_pair);
hold on;
plot(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',[plotcol 's'],'MarkerFaceColor',plotcol);
set(gca,'FontSize',16);
grid on
%hold on;
% %axis([0 0.5 0 200]);
h=xlabel('A/m');
% % %h=xlabel(['Longitude Bin ' sprintf('%4.1f ',ki) '^\circ  '  sprintf('%4.1f ',ki+incrlat) '^\circ']);
set(h,'FontSize',16);
ylabel('\DeltaH');
[c,sfit] = polyfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)',1);
m = polyval(c,[min(IndiaE(:,eej_index)),max(IndiaE(:,eej_index))]);
% %
% %plot(obst,chpt,'r.');
% % % xcor1(nloop-1) =dummy(2,1);
% [aaa,rstat] = robustfit(IndiaE(:,eej_index),obsE(:,obs_index).*(f_ch./f_ob)');% with added Zero the fitted curve is made to pass through origin
% dummy = corrcoef(rstat.resid,(obst-chpt)*24);
plot([min(IndiaE(:,eej_index)),max(IndiaE(:,eej_index))],m,plotopt,'LineWidth', 2);
text(0,25,sprintf('dH = %5.1f * I + %5.2f',c(1),c(2)),'FontSize',16,'Color',plotcol);
% end;
% text(0,20,sprintf(' CC %4.2f',cof(2,1)),'FontSize',16);
% %text(0,110,sprintf('dH = %5.1f * I + %5.2f',c(1),c(2)),'FontSize',16,'Color',plotcol);
% h=title(sprintf('%d ',x1data(nbin)));
% set(h,'FontSize',16);
% %end;
title(obs_title(i_pair,:));
axis([0 0.3 0 200]);
end;


%-----------------END
else,
    x1data(nbin) = NaN;
     
end;
nbin = nbin + 1;
clear ind f_ch f_ob ind1 ind2 ndata nnodata IndiaE obsE ltt1 ltt2 ltt3 ltt4 obst chpt X Y;
end;

x1data(isnan(x1data)) = [];
ndd(length(x1data)+1:end) = [];
ss= num2str(ndd');
fprintf('\n');
figure(1);
subplot(3,2,i_pair);
hold on;
errorbar(x1data,ceff,cerr/2,plotopt,'LineWidth',2);

text(x1data-2,ones([1,length(x1data)])+0.05,ss);
%axis([-80 80 -.6 1.4]);
grid on;
hold on;
text(x1data-2,ones([1,length(x1data)])+0.05,ss);
end;
more on